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THE  DETERMINATION  OF  CRITICAL  FLUTTER  CONDITIONS  OF  NONLINEAR  SYSTEMS 

by 

D.  L.  Woodcock 

SUMMARY 

A  flutter  analysis  procedure  for  nonlinear  systems  is  proposed  as  an 
alternative  to  timewise  integration  methods.  It  is  based  on  an  energy  method  due 
to  J.  Roorda  and  S.  Nemat-Nasser;  and  shows  promise  of  being  a  practical  pro¬ 
cedure  provided  the  number  of  parameters  can  be  minimised  by  the  representation 
of  the  flutter  system  by  an  equivalent  (condensed)  two  degree  of  freedom  system 
in  the  neighbourhood  of  the  critical  condition.  The  relationship  to  the  simpler 
method  put  forward  by  R.F.  Taylor  et  alia  is  considered. 
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1  INTRODUCTION 

The  theoretical  prediction  of  flutter  characteristics  can  normally  be  treated 
adequately  as  a  linear,  small  perturbation,  problem.  However  there  may  be  circumstances 
where  a  nonlinear  analysis  is  essential.  This  may  sometimes  be  the  case,  for  example, 
in  transonic  flow  conditions.  One  obvious  method  that  could  be  used  would  be  to  determine 
by  numerical,  timewise,  integration  the  motion  consequent  upon  a  number  of  initial  dis¬ 
turbances.  So  by  interpolation  one  should  be  able  to  determine  possible  states  of 
periodic  motion.  Alternatively  one  can  attempt  to  determine  directly  such  limit  cycles. 
This  present  paper  is,  therefore,  an  initial  consideration  of  how  this  might  be  done. 

An  approximate  method,  with  two  possible  simplifications,  is  suggested. 

2  THE  BASIC  RELATIONSHIPS 

Suppose  we  have  an  aeroelastic  system  where  Lagrange  equation  of  motion  can  be 
approximated  sufficiently  well,  for  the  purpose  of  finding  periodic  solutions,  by  the 
nonlinear  matrix  differential  equation:- 

Av2q"  +  ^Bv  +  ^  q'  +  ^C  +  q 

+  f(q,vq’)  +-^g(q,vvq’)  -  0  .  (2-1) 

v 

Here  A  is  the  inertia  matrix,  B  and  D  the  linear  aerodynamic  and  structural  damping 
matrices  respectively,  C  and  Ij  the  linear  aerodynamic  and  structural  stiffness 
matrices  respectively,  f  is  a  column  vector  of  the  nonlinear  aerodynamic  terms,  and 
g  is  a  column  vector  of  the  nonlinear  structural  terms.  A  non-dimensional  time  t  (=  u>t) 
has  been  introduced  where  oj  is  an,  as  yet,  undetermined  frequency.  The  frequency 
parameter  based  on  this  frequency,  the  airspeed  and  a  chosen  reference  length  is 
denoted  by  v  .  The  primes  denote  differentiation  with  respect  to  t  .  v  is  an 
airspeed  parameter,  being  the  ratio  of  the  airspeed  to  some  reference  speed.  The 
matrices  A  and  E  will  both  be  symmetric  and  positive  definite,  but  the  other  square 
matrices  (B,  C  and  D)  have  no  special  properties.  When  desired  we  will  divide  them  into 
symnetric  and  skew-symmetric  parts  denoted  respectively  by  the  subscripts  s  and  a  . 

We  wish  to  find  a  relationship  between  the  amplitude  of  steady  state  oscillations 
and  the  speed  parameter  v,  and  also  to  determine  whether  the  steady  state  oscillations 
are  stable.  Of  course,  one  may  wish  to  vary  another  parameter  rather  than  v  but  the 
procedure  should  be  basically  similar  to  that  which  we  will  describe.  It  is  not  to  be 
expected  that  steady  state  oscillations  will  be  possible  at  all  values  of  v  .  The 

method  to  be  described  is  an  application  of  that  proposed  by  Roorda  and  Nemat-Nasser ' . 

2 

The  rather  simpler  approach  suggested  by  Taylor,  Bogner  and  Stanley  ,  which  is  a 
generalisation  to  a  multi-degree  of  freedom  system  of  the  text  book  energy-balance  method 
(see  eg  Ref  3,  p  100),  will  appear  as  a  by-product. 

Let  us  assume  the  equation  (2-1)  has  a  periodic  solution  and  let  the  parameter  w 
which  we  have  introduced,  but  not  specified,  be  the  frequency  of  this  motion.  We  will, 
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as  an  approximation  to  the  solution,  take 

q  -  ap(x,T)  (2-2) 

where  a  is  an  amplitude  parameter,  and  x  is  a  column  vector  of  parameters  (i  =  l-*-m) 

which  are  to  be  determined  so  that  (2-2)  is  close  to  the  true  solution,  p  will  thus  be 
periodic  of  period  2ir  and  in  particular 

p(x,0)  =  p(x,2tt)  .  (2-3) 

The  existence  of  a  steady  state  solution  implies  that  no  energy  is  accumulated  or 
dissipated  over  a  complete  cycle.  This  means,  writing  the  left-hand  side  of  (2-1)  as 
y(q.q'.q")  »  that  our  assumed  solution  must  satisfy 

2ir 

a  I  |p'  (x,T)}Ty(ap,ctp'  ,ap")dT  =  0  .  (2-4) 

0 


The  matrices  A  ,  E  and  C  since  they  are  symmetric,  and  the  matrices  B  and  D 

s  a  a 

since  they  are  skew-symmetric,  will  make  no  contribution  to  this  integral.  Equation  (2-4) 
can  therefore  be  rewritten 


x(o,v,v,x) 


ap'  +  C  ap  +  f (ap.vap')  +  — ^  g(ap,vvap' )[dr  - 


0 

(2-5) 


where  we  have  assumed  a  is  not  zero. 


To  obtain  further  useful  relationships  we  turn  to  Hamilton's  principle  which 
states  that  the  victual  work  done  by  the  generalised  forces  during  any  admissible  virtual 
displacements  over  an  arbitrary  period  of  time  must  be  zero.  Now  the  elements  of  y  are 
the  generalised  forces,  since  (2-1)  is  the  equation  of  those  forces  to  zero,  and  so 


/ 


T 


1 


<5q  y(ap.ap,,ap")dT 


0  . 


(2-6) 


But  from  (2-2)  we  can  write  6q  in  terms  of  arbitrary  variations  of  the  elements  of  x 
and  of  the  frequency  w  .  Thus,  remembering  that  dt/dw  -  t  «  t/w  , 


6q 


p'  (x,t)6oj  + 


(2-7) 
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Consequently,  since  5w  and  Sx  are  arbitrary,  and  taking  the  time  interval  (tj  -*■  t 2) 
to  be  a  complete  cycle,  we  have  the  set  of  equations 


2tt 


T 

(x,x) |  y(ap,ap' ,ap")TdT  =  0 


2ir  \T 

/  pFrj  y(“p»“p’»«p")dT 


i  =  1  -*•  m  . 


We  have  here  assumed  that  a/w  ,  as  well  as  a  ,  is  not  zero. 

The  above  equations  (2-5),  (2-8)  and  (2-9)  provide  a  means  of  determining  all  but 

one  of  the  unknowns  a,  v,  v,  x  in  terms  of  the  remaining  one  -  say  determine  v,v 

and  the  x^  in  terms  of  the  amplitude  parameter  a  .  The  accuracy  of  the  approximate 

solution  will  depend  on  how  adequately  the  assumed  form  for  the  mode  of  displacement 

2 

p(x,t)  can  represent  the  true  solution.  The  method  of  imposed  disturbances  (or  energy- 
balance  method  (Ref  3,  p  100))  uses  just  equation  (2-5);  the  only  unknown  allowed  in  the 
assumed  displacement  is  the  amplitude  parameter  a  ,  the  frequency  and  displacement  mode 
being  determined  from  the  limiting  (linear)  case  when  a  0  .  An  alternative  simplifi¬ 
cation  would  be  to  just  use  equations  (2-5)  and  (2-8)  and  so  determine  say  the  frequency 
and  speed  parameters  in  terms  of  the  amplitude  parameter.  However,  as  the  example  of 
section  2.1  shows,  it  is  questionable  whether  this  would  yield  any  improvement  over  just 
using  (2-5) ,  though  it  does  avoid  the  assumption  that  the  frequency  is  unchanged  from 
the  linear  case. 

Before  further  consideration  of  how  this  suggested  method  could  be  implemented  we 
will  in  the  next  section  consider  a  very  simple  example  in  order  to  get  a  feel  for  what 
is  involved.  ) 

2.1  A  simple,  one  degree  of  freedom,  example 


Consider  the  equation 


v  x  +  —  (x2  -  v)x’  +  — =■ 
v  2 

v 


(2-10) 


This  represents  a  one  degree  of  freedom  system  in  which  the  structural  damping  is  the 
nonlinear  term  evx2x'/v  ,  and  the  only  aerodynamic  term  is  the  damping  (-evx')  .  By 
a  change  of  variables  this  can  be  transformed  into  the  standard  van  der  Pol  equation 


±4  *  uCE2  -  »  #  *  S 


(2-11) 


and  so  we  see  (of  Ref  3,  pp  102-106)  that,  for  e  small,  equation  (2-10)  has  a  limit 
cycle  solution 


6 


2/v  cos  t  +  0(e) 


with 


1  2 
v  -  -  +  0(e) 
v 


(2-12) 

(2-13) 


Note  that  the  solution  of  van  der  Pol's  equation  transformed  to  our  variables  gives 

x  =  2/v  cos | -  t|+  0(e)  ,  but  our  choice  of  w  as  the  frequency  of  the  periodic 

solution  means  that  the  period  in  t  is  2w  and  so  we  obtain  (2-13). 

Now  if  we  neglect  the  nonlinear  term  in  (2-10)  we  find  that  the  only  periodic 


solution  occurs  when  v 


and  vv  is  finite  and  is  (a  is  arbitrary) 


=  a  cost— I  . 
\vv/ 


(2-14) 


With  the  proviso  about  w  this  becomes 


with 


X  -  a  cos  t 


(2-15) 

(2-16) 


We  will  therefore  use  equation  (2-15)  as  the  trial  solution,  corresponding  to 
equation  (2-2),  for  the  methods  proposed  in  section  2  when  either  just  equation  (2-5)  or 
equations  (2-5)  and  (2-8)  are  used.  Equation  (2-5)  is 


which  gives 


2tt 
aev  f 

v  i 


.  2  2 
(a  cos  t 


v) 


.  2 
sin  t 


dr 


v 


In  addition  for  this  system  equation  (2-8)  becomes 


(2-17) 


(2-18) 


2ir 

-/I- 


2  ev  2  2  .  cos  t  ■  .  . 

v  cos  t  — —  (a  cos  t  -  v)  sm  t  +  - ^ —  r  T  smT  dt 


(2-19) 


and  so 


ip 


0  . 


(2-20) 


The  solution  of  this  equation  is 


v 


(-£♦  >) 


(2-21) 
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which  combined  with  (2-18)  gives 

v  -  I  +  0(e2)  .  (2-22) 

Thus,  in  either  case  the  solution  thus  obtained  has  an  error  0(c)  in  x  and  one  of 
2 

0(e  )  in  v  {of  equations  (2-1 2)  and  (2-13);  for,  when  only  (2-5)  used,  v  is  assumed 
to  be  given  by  the  linear  solution  (2-16). 

Making  use  also  of  (2-8)  avoids  the  assumption  that  the  frequency  is  unchanged  from 
the  linear  case,  but  otherwise  produces  no  improvement  and  it  is  obvious  that  this  must 
be  a  consequence  of  the  trial  solution  (2-15)  being  too  crude. 

To  obtain  a  store  accurate  solution  we  take  as  a  trial  solution  equation  (2-15) 
with  the  addition  of  a  third  harmonic,  ie 


x  «  a|cos  t  +  Xj  cos  3t  +  sin  3t|  =  ap  . 


(2-23) 


This  yields,  omitting  second  order  terms  in  Xj  and 


and 


p'  -  -  sin  r  -  3Xj  sin  3t  +  3x^  cos  3t 


ev  .  ,  2  2  . 

cos  r  -  —  smr  (ct  cos  t  -  v) 


(2-24) 


y(ap,ap'  ,ap")  *  a  ^  j  ^  -  v2^ 

+  xij(7“ 9v2) C08  3t  ~t( 

*  *2^  ■  9v2)  8in  3t  ~  T  (2  si 


.  2  2 
2  sin  t  cos  t  cos  3t  +  3(a  cos  t  -  v)  sin 


2  2 

sin  x  cos  t  sin  3t  -  3(a  cos  t  -  v)  cos  3t] 


i 


,(2-25) 


Thus  we  obtain  from  equations  (2-5),  (2-8)  and  (2-9)  to  the  same  order  of  accuracy 


I  jsin  t  (a  cos  T  -  v)  +  xj|2  sin2r  cos  t  cos  3t  +  6  sin  t  sin  3t  (a2  cos2t  -  v)^ 

0  L 

4  x2^2  sin2t  cos  t  sin  3t  -  6  sin  t  cos  3t  (a2  cos2r  -  v^Jdt 


-  0  . 
(2-26) 
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2* 


[  (?  ’  ”1  “ 


.  e\>  .  2  2  2 

sin  t  cos  t  — —  sin  r  (a  cos  r  -  v) 


+  x. 


•'1(7-') 

ev  /.  . 

- (  2  si 

v  \ 

(V’”2)51 

-v( 


sin  r  cos  3t  +  3| 


ft  ■  ■') 


cos  t  sin  3t 


2  2  2 
sin  t  cos  t  cos  3t  +  6  sin  t  sin3t  (a  cos  t 


-)! 


sin  t  sin  3t  -  3 


ft  -  •') 


cos  T  cos  3t 


m  2  2  2 
2  sin  t  cos  t  sin  3t  -  6  sin  t  cos  3t  (a  cos  t 


-)( 


TdT  =  0  .  (2-27) 


2tt 


l  0-1 


G  V  2  2 

cos  t  cos  3x  -  —  sin  t  cos  3t  (a  cos  x  -  v) 


+  x 


lft-1 


2 

cos  3t 


ev  /„  2  22 

"  ~  \2  Sln  T  cos  T  cos  3t  +  3(a  cos  t  -  v)  sin  3t  cos  3t 


)| 


! 
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Performing  the  integrations  we  get  respectively 


(4 -’)**■(-  7*  ¥)  ■  0 

1  /  I  2\  nev  /a2  \  \  2  ucv  foa2  -  l\(  .  (  ev  ( 7  a2  3v\! 

■ 2 17 - v )  - ~  b-  ■  vJ  *  *>  r  7  —  IrHl  l2r~  fa*T-r)| 


(4  -  9v2)xi  +  ^  (t  - v)  x2  -  0 

eva2  3 ev  /V  \  /  I  2\ 

TT  •  v  [  2  “  V)  *1  l  2  "  9V  ) 


x2  =  0 


The  solution  of  these  equations  is 


3a4  2  nf  3y 

-53-  e  +  0(e  > 


x2  =  -  §2  6  +  0(e3) 


a2  .  3a4(3a2  -  1)  2  3V 

V  =  T+ - i"28 -  £  +0(E  ) 


4  a2  (.  la2  a4\  2  .  3. 

2  "  8  \3  768  1%)  6  0(e  5 


and  so  from  (2-23) 


3 

x  ■  a  cos  t  — sin  3x  +  O(e^)  .  (2-38 

.  .  4 

Comparing  with  the  known  solution  ( eg  McLachlan  ,  p  45)  it  is  easily  seen*  that  we  have 

,  2 
now  got  a  solution  whose  error  in  both  x  and  v  is  of  0(e  )  . 

3  THE  METHOD  OF  IMPOSED  DISTURBANCES  (ALIAS  ENERGY  BALANCE) 

When  only  equation  (2-5)  is  used  we  take  in  the  trial  solution  (2-2) 


p(x,r)  -  5  cos  t  -  n  sin  t 


*  Probably  the  easiest  w  •»  to  do  thi  is  to  put  v  ■  1  and  then,  from  (2-36), 

2  /  33  2\  .2 

0  "  M*  "  T  e  /  ’  --*• * 1  Cr  c'ie  coefficient  of  e  in  the  expression  for  v 

(equation  (2-37)) It  (9S3/4o#)  compered  with  the  correct  value  of  (-1/16). 
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where  (£  +  in)  is  the  right  hand  eigenvector  satisfying 


.  2 
Av*  + 


(3-2) 


in  the  critical  state  when  and  are  real.  Thus  v  ,  and  (£  +  in)  are 

obtained  by  the  solution  of  the  normal  linear  flutter  problem.  Two  elements,  all  told, 
in  5  and  n  will  be  arbitrary.  We  also  assume  that  the  frequency  in  the  nonlinear 
case  is  the  same  as  in  the  linear  case,  ie  we  take 


vv 


(3-3) 


Other  assumptions  could  be  made  but  in  the  absence  of  any  strong  evidence  the  above 
seems  a  reasonable  choice.  Substituting  from  (3-1) and  (3-3)  in  (2-5)  and  integrating 
the  linear  terms,  then  gives  the  following  equation  relating  the  amplitude  a  of 
maintained  oscillations  with  the  speed  parameter  v  . 


air<v 


*7  ■>.)*♦  "’(i*.  *7°,)" 


*  «  V 


2tt 

•  / 


T  T 

(£  sin  t  +  n  cos  t) 


otvrtv„ 

Z  Z 

a(5  cos  t  -  n  sin  t),  -  (-  5  sin  t  -  n  cos  t) 


—  g Jct( C  cos  T  -  n  sin  t),  av^v^I-  £  sin  T  -  n  cost)! 

. (3-4) 


dx. 


Alternatively  it  could  have  been  written  as  an  equation  connecting  a  and  the  frequency 
parameter  v  .  If  the  nonlinear  aerodynamic  and  structural  terms,  f  and  g  ,  can  be 
evaluated  numerically  with  little  difficulty  -  one  might,  for  example,  have  an  analytical 
approximation  to  g  deduced  from  resonance  tests  -  then  the  solution  of  (3-4)  should  be 
no  problem.  However,  one  may  obtain  f  or  g  by  a  timewise  integration  procedure. 

Thus  one  would  have  to  assume  values  of  a  and  v  (or  equivalently  v) ,  evaluate  the 
right  hand  side  of  (3-4)  and  then  compare  it  with  the  left,  and  then  repeat  for  other 
(a,v)  until  one  had  got  equality.  In  such  a  case  a  scalar  parameter  could  be  put  on, 
say,  Dg  (the  symmetric  part  of  the  structural  damping  matrix)  and  the  results  inter¬ 
preted  as  values  of  this  parameter  as  a  function  of  a  and  v  .  The  curve  in  the  (a,v) 
plane  where  the  parameter  was  unity  would  then  indicate  the  critical  states  -  the 
limit  cycles. 


3. I  Energy  balance  with  frequency  determination 

In  this  section  we  consider  in  more  detail  the  use  of  (2-8)  along  with  (2-5).  In 
the  trial  solution  p(x,x)  is  again  taken  to  be  given  by  equations  (3-1)  and  (3-2).  The 
assumption  (3-3)  is,  however,  no  longer  needed.  Thus  (2-5)  gives 


airlv? 


+  2^ 


(b  +  -  D  +  viiT(b  +  -  d  )n 

\  S  V  8/  \  S  V  S  / 

2tt 

=  j  sin  t  +  nT  cos  t)  |  f  [a ( C  cos  x  -  n  sin  r),  ctv(-£  sin  T  -  n  cos  x)J  +  — j  x 

0 

xg|a(C  cos  x  -  n  sin  x),avv(~5  sin  x  -  n  cos  x)jjdx  (3-5) 


while  from  (2-8)  we  have 


f 

1  TJ 

f  2 

1 

a] 

2tt  C  -  irv 

B  +  -  D 

1  ' 

^  3 

s  v  sj 

r  i 

TT 

1  „  2.’ 

B  +  -  D 

+  -T 

C  +  — E  -  v  A 

l_s  vs. 

2 

L s  v2  J 

(b  +  -  D 

-  4 

C  +  E  -  v2a! 

Ls  v  s_ 

2 

L  s  v2  J 

=  j  (£T  sin  x  +  nT  cos  x)|f  Ja(C  cos  x  -  n  sin  x),  otv(-£  sin  x  -  n  cos  x)J 
0 

+  ~  g|a(?  cos  x  -  n  sin  x) ,  avv(-£;  sin  x  -  n  cos  x)J 


xdx 


,(3-6) 


These  two  equations  are  particularly  convenient  in  the  case  when  the  linear  structural 
damping  (or  indeed  just  the  symmetrical  part  of  the  linear  structural  damping  matrix) 
is  zero  and  there  is  no  structural  nonlinearity.  It  is  common  practice  to  assume  no 
structural  damping  and  rely  on  what  there  is  present  in  practice  for  an  extra  little  bit 
of  safety  margin*.  In  this  particular  case  Dg  and  g  will  be  zero  therefore.  The 
first  equation  (3-5),  which  then  becomes 

air|v|cXBg5  +  nTBgn)  +  24TCan 

2xt 

(£X  sin  x  +  nT  cos  x)f[a(5  cos  x  -  n  sin  x),otv(-£  sin  t  -  n  cos  x)J  dx  (3-7) 

0 

can  be  solved,  on  the  lines  discussed  in  section  3,  for  v  as  a  function  of  a  . 


*  The  addition  of  structural  damping  does  not,  however,  necessarily  make  a  system  more 
stable. 
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The  speed  parameter  v  is  then  immediate 1  •  evaluated  from  equation  (3-6)  which  in  this 
particular  case  gives 


2  T  T 

v  (C  EC  -  n  En) 


CT(v2A  -  2itvBg  -  Cg)c  -  nT^v2A.  +  2itvBg  -  Cgjn 

2ir 

+  —  J  (CT  sin  t  +  cos  x)f  [a(C  cos  x  -  n  sin  x) , 

0 

av(-C  sin  t  -  n  cos  x)Jxdxl  .  (3-8) 


When  the  only  nonlinearity  is  structural  the  solution  is  in  one  sense  rather  more 
difficult.  We  can  eliminate  v  between  the  two  equations,  (3-5)  and  (3-6),  and  obtain 
an  equation  connecting  a  and  w  .  Thus  writing  the  following  four  constants  as 


kl 

T 

=  C  B  C  + 
s 

T 

n  b  n 
s 

(3-9) 

T 

T  1  T 

(3-10) 

k2 

=  C  B  C  + 
s 

n  b  n - B  n 

S  TT  S 

C1 

T 

=  2CTCan 

(3-11) 

C2 

=  2CTC  n 
a 

*  -  "V} 

(3-12) 

and  also  using  the  notation 


g,  (ct.vv) 


g2(a,vv) 


- 2~  /  (CT  sin  t  +  nT  cos  x)g[ct(C  cos  x  -  n  sin  x) , 

ot(vv)  tt  J 

avv(-  C  sin  x  -  n  cos  x  )]  dx 

-  — (cXD  C  +  nTD  n) 

vv  \  s  s  / 

2tt 

- 1  v-  /  (CTsinx  +  nT  cos  x ) g  [a ( C  cos  x  -  n  sin  x), 

oc(vv)  *  * 

avv(-  C  sin  x  -  n  cos  x)Jxdx 


I 


(3-13) 


2ir(vv)‘ 


(CTEC  -  nTEn)  -  ±  (cTDgC  +  nTDgn  -  i  cTDsn)  +  ±  (cTac  -  nTAn) 


,(3-14) 


we  obtain  the  relationship 

(clk2  "  c2kl){k28l(a,vv)  “  k|82(a,vv^  “  {c2gl(a*vv)  “  c182(a‘vv)}2  .  (3-15) 
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Having  solved  this  equation  for  vv  as  a  function  of  a  one  can  then  return  to 
equation  (3-5)  or  (3-6)  to  determine  v  and  hence  v  .  In  this  case  (no  aerodynamic 
nonlinearity),  with  the  notation  introduced  in  equations  (3-9)  to  (3-14),  they  become 


2 

v  gj  (a»w)  -  vkj  -  Cj  =  0  (3-16) 

2 

v  g2(a,vv)  -  vk2  -  c2  -  0  .  (3-17) 


4  THE  FULL  PROPOSAL 

Our  full  proposal  involving  the  use  of  equations  (2-5),  (2-8)  and  (2-9)  permits 
the  introduction  of  some  unknown  parameters  xr  -  the  elements  of  x  -  in  the  mode 
p(x,t)  of  the  trial  solution  (2-2).  In  the  one  degree  of  freedom  example  considered  in 
section  (2-1)  these  were  taken  to  be  coefficients  of  harmonics  higher  than  the  funda¬ 
mental  (of  equation  (2-23).  For  the  multi-degree  of  freedom  system  there  is  more  choice. 
One  could  alternatively  allow  variation  in  the  relative  values  of  the  elements  of 
p(x,t)  rather  than  introduce  higher  harmonics,  or  perhaps  use  some  mixture  of  the  two 
ideas. 

Baldock^’^  has  demonstrated  how  one  can  often  find  an  equivalent  two  degree  of 
freedom  system  which  will  represent  the  behaviour  of  a  system  with  many  more  degrees  of 
freedom  with  good  accuracy  at  or  in  the  vicinity  of  the  flutter  condition.  Let  the 
modes  of  these  two  degrees  of  freedom  be  represented  in  terms  of  the  original 
generalised  coordinates  by  the  column  vectors  and  6  ,  and  let  the  flutter  mode  of 
the  linear  system  be 


vO 


q  m  (zj  +  i)<f>  +  0  +  iz2)e  . 


(4-1) 


Then,  on  the  basis  of  this,  it  seems  reasonable  to  take  as  a  trial  solution  for  the  non¬ 
linear  system  , 

q  =  ap(x,T)  =  a[{(zl  +  *!>♦  +  e}  C09  *  -  ♦  (z2  +  x2)6}  sin  rj  .  (4-2) 

Substituting  this  expression  in  equations (2-5)  and  (2-8)  gives  the  equations  given 
in  the  previous  section  (3-1)  -  equations  (3-5)  and  (3-6)  in  the  general  case  -  with 


(z  j  +  x  j )  4>  +  e 


♦  +  (z2  +  x2)e 


(4-3) 

(4-4) 


Of  course,  we  only  wish  to  retain  first  order  terms  in  x.  and  x,  ,  and  so  writing 


«> 


(n) 


Xj-0 


x2-0 


(4-5) 

(4-6) 


f0(a>v,T)  -  f[o(C0  cos  t  -  n0  sin  r),av(-  sin  t  -  nQ  cos  t)J 
g0(“»^v,T)  =  gJa(5Q  cos  t  -  nQ  sin  r),avv(-^Q  sin  r  -  n0  cos  t)] 


equation  (3-5)  becomes 


"'oK  *  "  D-)£o  *  V"5K  *  i  \)\  *  *  2-|  [»«JK  *  i  °>  -  ^c.*] 

*  ^[’’"SK  *  i »,)» *  «k9]  | 

2ir 

=  j  (Cq  sin  t  +  rip  cos  x)|f0(a,v,T)  4  -L  go(a,v  v,T)|dT 


ZTT 

+  X1  /  [<t’T  sin  T  f0(“*^T)  +  (c0  sin' T  +  no  cos  T)  * 

0 

*  {FJq> («.>». T)*a  cos  T  -  F^q  ^(o.v.t)^  sin  tjJdT 
27 r 

+  x2  /  [f  cos  T  f0(a*v,T)  ~  (e0  s^n  T  +  no  cos  T)  * 

o 

*  {^q)(a,v,T)6(>  sin  t  +  F^q  ^(a,v,r)0a  cos  rjjdT 
2tt 

+  ~  j  J^1  sxn  x  Sg^a,vv,T^  +  sin  t  +  Ti^  cos  t]  x 

x  |Gq^^  (a,vv,T)i>a  cos  t  -  G^q' 5  (a.vv.THa  sin  T|JdT 

x  2tt  f 

+  ~  J  [_9  cos  T  g0(a.vv,r)  -  sin  r  4  cos  rj  * 

X  {G^q)(a,vv,T)ea  sin  T  +  G^  ^(a,vvfT)0a  cos  tJ>JdT  (4 


where  {of  section  4.1)  \  G^q  ^  are  square  matrices  whose  jth  columns 

are  respectively 


S© 

00 

O' 


’j  L 


1  3qj 


q=a(£0  cos  x-r>0  sin  t) 


af(q,  Jq')l 

J  J  q=a(£g  cos  T-nQ  sin  t) 


V 


9gfq.»  wq' )  I 

3qi  I 

J  Jq-a(5Q  cos  T-n0  sin  t) 


and 


sq-  i 

J  Jq«a(SQ  cos  T-nQ  sin  t) 


(4-10) 


q^  being  Che  jth  element  of  q  .  We  are  assuming  that  the  functions  f  and  g  have 


a  Taylor  expansion  at  this  point.  Similarly  equation  (3-6)  becomes 


-  irv 


E  -  v  A 


[B.  *  i  *  £?(*2»[ss  *  i D.]  *  i[c.  *  7 

*v»J  -  f  [c. " 7 E '  ''2a])"o 


+  x. 


2”2'lC  a  +  nl{ 

^-4 

„  2  T  T/ 

1  _  2  A" 

2*  *0C.  +  ">0( 

Cs+7e-va) 

2ir 


=  j  sin  t  +  cos  T){f0(a»v>T)  +  g0(ct,vv,T)|tdT 


2»  r  t 

I  U  sin  t  fQ(a,v,T)  + 

(?q  sin  t  +  rig  cos  xj 

|  (a,v,t)  $a  cos  T  -  F^q  \a,v,T>4>a  sin  t|  rdx 


2tt 


x2  J  |eT  cos  t  fQ(a,v,T)  -  ^  sin  t  +  rig  cos  tJ  x 

x  {F£q)(«,v,T)ea  sin  r  +  F ^  ^(a,v,T)6a  cos  T^rdT 


2n 


*  ”5  /  L  sin  T  80(a*vv*T)  +  K  8in  T  +  n0  cos  T)  * 

x  {c£q)  (a,vv,t)*o  cos  t  -  G^q  ^(a,vv,T)$a  sin  tlltd-c 


0 


+  ~  J  L  cos  T  8o^ot,vv,T^  “  ro  sin  T  +  no  cos  T)  * 

x  (a,vv,r)0a  sin  t  +  G^q  Ha,vv,x)0a  cos  xjfj Tdr  . 


(4-11) 


Furthermore,  from  equation  (2-9),  with  p(x,t)  given  by  (4-2),  we  obtain  the  following 
two  equations,  keeping  again  only  first  order  terms  in  and  Xj  : 


am<(> 


r< 

1  /  E  2\ 

/  D\  | 

1 

/  E  2\ 

/ 

„  D\  a 

j 

i 

lr7“u)5»' 

V(B  *vj  n0| 
/ 

>  +  X, 

1  1 

(cs +  7  ‘ Av ) 

♦  -  vx2  ( 

B  +  -  6 

vi 

2ir 


cos  T 


+  ax. 


=  -  4>r  J  |{fo(a»v»T)  +  ~7  g0(“.yv.T)} 

«,{FJ‘l)(afv,T)*  cos2x  -  Fq*1  ^(a,v,xH  sin  t  cos  xj^ 

— (a,vv,x)$  cos2t  -  Ggq  Ha,vv,x)<t>  sin  t  cos  t| 
v  ' 

ax2|F^q^ (a,v,x)e  sin  t  cos  t  +  F^q  Ha,v,x)0  cos2t| 

— r{Goq)(a«VV’T)0  S^n  T  C0S  T  +  Goq  Ha,vv,x)0  cos2x[ 

V  ' 


dx 


.(4-12) 


and 


ax0T  He  +  ~  -  Av2 


2w 


=  0T 


!)n0  +  V(B  +  ?)  Co|  +  VXl(B  +  £)  *  +  Xl(Cs  +  4  “  Av2) 

I  [jfo(a’v’T)  +  ~T  8o(a,VV,T)}  sin  T 

j^F^q^ (a,v,x)$  sin  x  cos  x  -  F^q  Ha,v,x)$  sin2x| 


+  ax, 


+  — j-|GQqHa,vv,x)$  sin  x  cos  x  -  G^*1  Ha,vv,x)4>  sin2xi 
v  '  f 

-  ax2|FQqHoi,v,x)0  sin2x  +  F^q  Ha,v,x)0  sin  x  cos  x| 

- ^•iG^qHa,vv,x)6  sin2x+  G^q  Ha,vv,x)6  sin  x  cos  xildx 


■  1 
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Thus  we  have  four  equations,  (4-9),  (4-11),  (4-12)  and  (4-13),  which  can  be  solved  for 
the  four  of  the  five  unknowns  a,  v,  v,  and  x 2  in  terms  of  the  other  one.  The 
equations  are  linear  in  the  modal  parameters  Xj  and  X2  but  not  in  the  other  unknowns. 

The  solution  of  these  equations  should  not  be  as  troublesome  as  it  might  appear. 
Consider  first  the  case  when  there  is  no  structural  damping  and  no  structural  non¬ 
linearity  (D  and  g(q.vvq')  are  both  zero).  Then  all  the  integrals  that  appear  are 
independent  of  v  and  so  we  can  assume  values  of  a  and  v  and  obtain  respectively, 
from  equations  (4-9),  (4-11),  (4-12)  and  (4-13),  equation  of  the  form 


y01Xl  +  y02X2  +  y03 


^W1 1  ~  ~^2  Y1 l^xl  +  (y12  ~  ^2  Y 1 2^x2  +  (u13  "  ^2  ^3) 

^21  ”  ^2  Y21^X1  +  y22X2  +  (y23  ~  "TT  ^3) 

(W32  “  ~1.  y3l)X2  +  (y33  "  ^2  Y33) 


=  0 


y31Xl  + 


0  . 


The  coefficients  will  all  be  functions  of  a  and  v  .  For  example 


y02  =  2a,r 


2-it 

[^n0Bs9  +  ^0Ca9]  “  j  |^T  cos  T  Va’V*T>  ”  (^0  sin  T  +  n0  C0S  T)  X 


(4-14) 


(4-15) 


(4-16) 


(4-17) 


^F^  (a,v,r)0a  sin  r  -  F^*5  ^(a,v,t)9a  cos  x|JdT  • 


,(4-18) 


The  other  coefficients,  the  ,  are  constants.  An  example  is 


y 


32 


-  Tt6TE0 


(4-19) 


The  three  equations  (4-15)  to  (4-17)  form  an  eigenvalue  problem  which  can  be  solved 

2 

by  the  usual  methods  for  a/v  ,  x.  and  x_  .  We  require  a  solution  for  which  all 

!  ^2 

these  quantities  are  real  and  the  first  one,  (a/v  ),  is  positive.  There  always  will  be 
at  least  one  real  solution.  The  values  of  Xj  and  x2  ,  from  suitable  solutions  of 
(4-15)  to  (4-17),  can  then  be  substituted  in  equation  (4-14)  to  see  how  closely  it  is 
satisfied.  By  repeating  the  procedure  for  other  assumed  values  of  a  and  v  one  should 
then  be  able  to  find  the  conditions  under  which  all  four  equations  are  satisfied;  thus 
obtaining  v,  x]f  x2,  and  v  as  functions  of  a  for  the  critical  (limit  cycle)  flutter 
condition. 

When  the  only  nonlinearity  is  structural  it  is  convenient  to  use  a  notation  which 
is  an  elaboration  of  that  used  in  that  section  3.1  (see  Table  1  for  full  details). 
Equations  (4-9),  (4-11),  (4-12)  and  (4-13)  then  become 
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2 

-  a  factor  air  or  an  has  been  deleted  as  convenient 

-  respectively 

v2{g10(ot,vv)  +  xl8ll(ol,vv)  +  x2gi2(a,vv)} 


w2{g20<a,vv)  +  xlg21 ^a*vv^  +  x2g22^a,vv^} 


v2^g30(oi,vv) 

+  x]g3)(a,vv) -x2g32(a,vv)| 

-v(k30  +  k32X2)  +  (c30  +  C31Xl)  * 

0 

(4-22) 

v2{g40("*vv) 

+  x,gA](a.w) -x2g42(a,vv)| 

‘V<k40  +  k41Xl)'Cc40  +  C42X2)  * 

0 

.  (4-23) 

v(k10  +  kllXl  +  kl2x2> 


(C10  “  C!IX1  +  c!2*2)  *  ° 


(4-20) 


v(k20  ♦  k2,x,  ♦  k22x2) 


(c20  ‘  C21X1  +  C22X2) 


(4-21) 


The  procedure  that  can  be  followed  will  then  be  similar  to  that  just  suggested  in  the  case 

of  aerodynamic  (and  not  structural)  nonlinearity.  Values  of  a  and  vv  assumed;  the 

2 

eigenvalue  value  problem  posed  by  three*  of  the  above  four  equations  is  solved  for  v  , 

Xj  and  x2  ;  the  solution  is  substituted  in  the  remaining  equation;  and  the  process 
is  repeated  with  other  assumed  values  of  a  and  vv  until  the  locus  of  critical  condi¬ 
tion  is  found. 

4.1  The  matrices  etc 

No  difficulty  should  be  encountered  in  evaluating  the  matrices  whose  elements  are 
defined  by  equations  (4-10)  if  analytical  expressions  are  known  for  the  column  vectors 
f(q.vq’)  and  g(q,vvq').  There  may,  however,  be  cases  when  all  that  is  available  are 
numerical  values  Of  the  vectors  at  instances  during  a  specified  motion.  For  example 
this  most  likely  will  be  the  case  for  f  as  given  by  nonlinear  transonic  flow  calcula¬ 
tions.  At  first  sight  one  may  imagine  in  such  a  case  that  we  are  left  with  an  almost 
overwhelming  problem. 


However,  let  us  first  note  that  what  is  required  in  equations (4-9)  and  (4-11)  to 
)  is  not  t 
column  vectors 


(4-13)  is  not  the  two  square  matrices  F^  and  ^  in  isolation  but  the  two 


$  cos  t  —  F<q,)*  sin  r)  *  ± 


q-ap 


and 


(f^Q  sin  t  *  F^'h  cos  t)  -  -i  (~\ 

/(5-»Pr 


(4-24) 


(4-25) 


*  Stability  considerations  (of  section  4.1)  suggest  that  it  may  be  most  instructive  to 
take  the  last  three  -  (4-21)  to  (4-23),  A  solution  which  gives  real  and  positive 
(and  hence  x{  and  x2  real)  is  required. 
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where  pq(t)  =  p(0,t)  .  (4-26) 

Thus  by  evaluating  f(q.vq')  for  q  =  ap^  ,  and  for  one  or  two  adjacent  q  caused  by 
small  variations  in  Zj  and  z^  ,  one  should  be  able  to  evaluate  the  required  vectors. 

4.2  Stability  of  limit  cycles 

When  the  function  x(“«v,v,x)  (equation  (2-5))  is  positive*,  it  means  that  energy 
has  to  be  supplied  to  the  system  to  maintain  the  assumed  motion.  Our  proposed  solution 
procedure  determines  a  curve,  in  say  the  (a,v)  plane,  where  X  =  0  ;  and  also  provides 

values  of  x  at  other  points  in  this  plane.  If  x  is  negative  above  this  critical 

curve,  ie  when  a  is  increased  by  6a  ,  or  positive  below  the  curve,  then  we  can 
certainly  say  that  the  limit  cycles  corresponding  to  points  on  the  curve  are  unstable 
for  they  are  unstable  with  respect  to  perturbations  in  a  .  We  cannot,  however,  say 
with  absolute  certainty  that  the  limit  cycles  are  stable  if  the  contrary  is  true, 

that  is  if  x  is  positive  above  and  negative  below  the  critical  curve.  This  is 

because  the  motion  may  be  unstable  with  respect  to  other  perturbations  (of  Ref  ))• 

However,  in  practice,  one  would  expect  the  signs  of  x  »  adjacent  to  the  critical  curve 
X  *  0  ,  to  be  a  good  enough  guide  to  the  stability  or  otherwise  of  the  limit  cycles. 

If  there  is  a  stable  limit  cycle  then  it  means  that  we  can  have  what  is  usually  known 
as  limited  amplitude  flutter;  while  if  there  is  just  an  unstable  limit  cycle  then 
catastrophic  (diverging  to  infinity)  flutter  will  be  possible  when  the  system  undergoes 
a  big  enough  disturbance. 

5  CONCLUDING  REMARKS 

The  question  as  to  whether  procedures  for  determining  limit  cycles,  such  as 
those  described  in  this  Memorandum,  are  more  attractive  than  other  approaches  will  depend 
largely  on  the  number  of  parameters  necessary  in  the  description  of  the  limit  cycle. 

If,  as  in  the  detailed  development  of  section  4,  one  can  get  away  with  two  parameters 

x.  and  x.  ,  in  addition  to  the  amplitude  and  frequency  parameters  (a  and  v)  (of  equation 

•  £  5  6 

(4-2)),  then  it  may  well  prove  the  more  economical  method.  Baldock's  way  *  of  con¬ 
densing  a  linear  flutter  system  to  an  equivalent  binary  certainly  provides  a  promising 
guide  to  the  choice  of  parameters  in  such  a  case. 


*  Or,  equivalently,  the  left-hand  side  is  greater  than  the  right-hand  side  in 
equations  (3-4)  or  (3-5)  or  (4-9),  or  the  left-hand  side  of  (4-20)  is  negative. 
The  energy  supplied  is  also  zero  in  the  trivial  case  a  ■  0  (of  equation  (2-4)). 
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Table  1 

NOTATION  USED  IN  EQUATIONS  (4-20)  TO  (4-23) 


;10 

2CSCan0  * 

C1 1 

2nOCa*  • 

'12  ■  2£SC.6  • 

20 

2£5c."o  *  i 

4£5Cs£0  -  "SC«"o}>  C21 

■  2"JC,»  -  7  • 

c  *  2^C  e--  ruC  0 

31  0  a  it  0  s 

30 

*?c£o  • 

C31  ’ 

T 

<P  CJ  , 

s 

40 

0TCno  , 

C42  " 

T 

e  c  e 

s 

k10  C0Bs50  +  n0Bsn0  * 


1  T 

1  k  A  n 


k30  =  *  Bn0  ’ 


k40  “  9  Be0  * 


:i  i 

-  2£5» 

:  i 

2CW 

T 

'32 

=  4>  B0 

T 

'41 

»  8  Bf 

ki2  2*sv  * 


k22  '  2'l*sd  ~  7  * 


g10<«.vv) 


g, , (o»VV) 


1 


air(vv)' 


2ir 

/  i^O  8in  T  +  n0  C0S  T)«0(a'VV*T)dT  "  (£0DS£0  +  nODsno) 
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constant  coefficients  in  equations  (4-20)  to  (4-23)  (no  aerodynamic 
nonlinearity,  of  Table  1) 

column  vector  of  nonlinear  aerodynamic  terms 
see  equation  (4-7) 
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